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We report a general class of steady and transient states of granular gases. We find that the 
kinetic theory of inelastic gases admits stationary solutions with a power-law velocity distribution, 
f(y) ~ v~ a . The exponent a is found analytically and depends on the spatial dimension, the degree 
of inelasticity, and the homogeneity degree of the collision rate. Driven steady-states, with the same 
power-law tail and a cut-off can be maintained by injecting energy at a large velocity scale, which 
then cascades to smaller velocities where it is dissipated. Associated with these steady-states are 
freely cooling time-dependent states for which the cut-off decreases and the velocity distribution is 
self-similar. 
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I. INTRODUCTION 

The statistical physics of granular gases is unusual in 
many ways 0, S 13] • Shaking a box of beads, no mat- 
ter how hard, fails to generate a thermal distribution 
of energy. Instead, the velocity distributions arc not 
Maxwellian 0, |M and energy may be distributed un- 
evenly in space 0,S or among different components of 
a polydisperse granular media |loL fTT| . Moreover, spatial 
correlations may spontaneously develop [lj. Granular 
gases also exhibit interesting collective phenomena such 
as shocks [ij, Q, clustering pi E [H E E3, and 
hydrodynamic instabilities |20l |21|. Energy dissipation, 
which results from inelastic collisions, is largely respon- 
sible for this rich phenomenology. 

Dilute granular matter can be studied systematically 
using kinetic theory. This approach has been used to 
quantitatively model situations where the dynamics are 
primarily collisional |22l l23l l2 ll |25| . Kinetic theory has 
been used to derive transport coefficients in the contin- 
uum theory of rapid granular flows, and it has also been 
used to model freely evolving and driven granular gases. 

Spatially homogeneous systems are a natural starting 
point for investigations of granular gases. Theoretical, 
computational, and experimental studies show that the 
system cools indefinitely without energy injection, and 
that it reaches a steady-state when energy is injected 
to counter the dissipation. In the freely cooling case, 
the velocity distribution follows a self-similar form and 
in the forced case, the velocity distribution approaches 
a steady-state. In either case, the velocity distributions 
have sharp tails, and in particular, all of their moments 
are finite. 

In this study, we consider the very same spatially ho- 
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mogeneous granular gases and show that there is an ad- 
ditional family of steady and transient states. First, 
we demonstrate that for a special, analytically solu- 
ble case, the unforced Boltzmann equation admits non- 
trivial stationary states where the velocity distribution 
has a power-law high-energy tail. Then, we show that in 
general, the tail of the distribution obeys a linear equa- 
tion and use this master equation to demonstrate that 
stationary states with power-law tails are generic, exist- 
ing for arbitrary dimension and arbitrary collision rules. 
The characteristic exponents are obtained analytically 




The mechanism responsible for these stationary states 
is an energy cascade from large velocity scales to small ve- 
locity scales that occurs due to the inelastic particle col- 
lisions. Driven steady-states with the same characteristic 
exponent and a high velocity cut-off can be maintained 
by injecting energy at a large velocity scale to compen- 
sate for the energy dissipated in the cascade. We confirm 
these steady-states using Monte Carlo simulations. We 
propose that such steady states can be experimentally 
realized in driven granular systems in which energy is 
injected at large velocities. 

There is also a family of closely related freely cooling 
time-dependent states. We demonstrate this explicitly 
in one-dimension. In these transient states, the velocity 
distribution coincides with the stationary distribution up 
to some large velocity scale, but falls off exponentially 
beyond that scale. This cut-off velocity obeys Haft's 
cooling law and decreases algebraically with time until 
the power-law range collapses. The velocity distribution 
is self-similar and the underlying scaling function is ob- 
tained analytically using the linear Boltzmann equation. 
These freely cooling states are confirmed using numerical 
integration of the Boltzmann equation. 

This paper is organized as follows. The system is set- 
up in section II and a special case is solved in section III. 
Dynamics of large velocities and the linear Boltzmann 
equation are described in section IV. Stationary states 
are detailed in section V, driven steady-states in section 
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VI and transient states in section VII. We conclude in 
section VIII. 



II. INELASTIC GASES 

We study a spatially homogeneous system of identical 
particles undergoing inelastic collisions. First, we con- 
sider one-dimension where the linear collision rule is 



"1,2 = pu\-i + qu 2 ,i 



(1) 



with v\ t 2 the post collision velocity and u\ t2 the pre- 
collision velocity. The collision parameters p and q obey 
p + q = 1 . The relative velocity is reduced by the resti- 
tution coefficient r = 1 — 2p as follows: (i>i — v 2 ) — 
—r(u\ — u 2 ). In each collision, momentum is conserved, 
but the total kinetic energy decreases. The energy loss 
is AE = pq{u\ — U2) 2 ■ Energy dissipation is maximal 
for the extreme case of completely inelastic collisions 
(r = Q,p = 1/2) and it vanishes for the extreme case 
of elastic collisions (r = l,p = 0). 

In this study, we consider the general collision rate 



K(ui,u 2 ) = \ui - u 2 \ 



(2) 



with < A < 1 the homogeneity index. For particles 
interacting via the central potential U(r ) ~ r~ u , the ho- 
mogeneity index is A = 1 — 2 [27I l2q . There are two 
limiting cases: (i) Hard-spheres, where the collision rate 
is linear in the velocity difference, A = 1, are used to 
model ordinary granular media; (ii) Maxwell-molecules, 
where the collision rate is independent of the velocity are 
used to model gra nular media with certain dipole inter- 
actions [HEMES!- 

Let f(v, t) be the distribution of particles with velocity 
v at time t. It is normalized to unity, J dvf(v) = 1 
(henceforth the dependence on t is left implicit). For 
freely evolving and spatially homogeneous systems the 
distribution obeys the Boltzmann equation 



dfjv) 
dt 



du!du 2 \ui - u 2 \ x f{ui)f(u 2 ) 



S(v — pui — qu 2 ) — S(v — U\) 



(3) 



In this master equation, the kernel equals the collision 
rate (J2J) and the gain and loss terms simply reflect the 
collision law QJ. The Boltzmann equation assumes per- 
fect mixing as the probability of finding two particles at 
the same position is taken as proportional to the prod- 
uct of the individual particle probabilities. It is exact 
when the strong condition of perfect mixing or "molec- 
ular chaos" is met, but it is only approximate when the 
particle positions are correlated. 

One well-known solution of the this equation is the 
"homogeneous cooling state" where the velocity distri- 
bution is self-similar in the long time limit [34l |35| , 



(4) 



with the characteristic velocity vq 
analysis, the collision rate K oc v, 
proportional to time, K ~ 
cooling law E 



Applying dimensional 
should be inversely 
This leads to Haff's 



V ~ t 



-1/X 



(5) 



Alternatively, it follows from the rate equation 
dvo/dt oc — i>o +A , implying that exponential decay occurs 
for the limiting case of Maxwell molecules. Statistics 
of energetic particles are characterized by the tail of 
the distribution, and for freely cooling states, there is 
a stretched exponential decay |3J, |3|| |23 



i/>(z) ~exp (-|z| A ) 



(6) 



for A > as \z\ — > 00. 

In the freely cooling states, all energy is dissipated from 
the system and the particles come to rest, f(v, t) — * S(v) 
as t — ► 00. Thus, the system reaches a trivial stationary 
state. Are there any nontrivial stationary states? Quite 
surprisingly, the answer is yes. Our main result is that 
generically, there is a family of nontrivial stationary so- 
lutions of the Boltzmann equation. 



III. AN EXACT SOLUTION 

The stationary velocity distribution can be obtained 
analytically for one-dimensional Maxwell molecules. 
Since the governing equation (J3J is in a convolution 
form, it is natural to employ the Fourier transform [39J, 
F(k) = J dve lkv f(v). The stationary state (d/dt = 0) 
satisfies the non-local and non-linear equation ^(]> S3 



F(k) = F(pk)F(qk). 



(7) 



Normalization implies F(0) = 1. 

For elastic collisions, p = 0, every distribution is a 
stationary state, but this is a one-dimensional anomaly, 
because in higher dimensions, the stationary distribution 
is always Maxwellian |29j. For all < p < 1 and p+q = 1, 
there is a family of stationary solutions 



F(fc) = exp(-« |/c|) > 



(8) 



characterized by the arbitrary typical velocity vq. Per- 
forming the inverse Fourier transform, the velocity dis- 
tribution is a Lorentz (Cauchy) distribution 



f(v) 



■kvq 1 + (v/v y 



(9) 



This distribution decays algebraically at large velocities. 
For freely cooling Maxwell-molecules in one-dimension, 
the velocity distribution has a related form, a squared 
Lorentzian [43|. 

This stationary distribution does not evolve under the 
collision dynamics since at each velocity there is perfect 
balance between collisional loss and collisional gain. The 
total energy density and the total dissipation rate are 
both divergent due to the shallow tail of the velocity dis- 
tribution. 



IV. CASCADE DYNAMICS 

To analyze the general behavior, we focus on the dy- 
namics of very energetic particles. This allows us to de- 
rive the power-law decay and to obtain the characteristic 
exponent for all spatial dimensions and all collision pa- 
rameters. 
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One-Dimension 



FIG. 1: The cascade process. 



The collision integral in Eq. Q greatly simplifies in 
the limit v — > oo. Since the distribution decays sharply 
at large velocities, the product f{u{]f(u2) is maximal 
when one of the pre-collision velocities is large and the 
other small. For the gain term there are two possibilities: 
either u\ 3> ui and then v — pu\ or ui 3> u\ and then 
v = qu 2 - Let us denote the large velocity by u and the 
small one by w. The double integral separates into two 
independent integrals, 



B. Arbitrary Dimension 



df(v) 
dt 



dwf{w) J du\u\ x f(u) 
x [5(v — pu) + S(v — qu) 



(10) 



5(u)]. 



Here, the collision rate \u 



w\ 



was approximated by 



\u\ . The integral over the smaller velocity equals one, 
and performing the integration over the larger velocity 
yields 



df(v) 
dt 



1 



P 



l+A 



/(- 

P 



(11) 



The tail of the velocity distribution satisfies a non-local 
but linear evolution equation. 

The linear Boltzmann equation is valid for broader 
conditions compared with the full nonlinear Boltzmann 
equation. The only requirement is that energetic par- 
ticles are uncorrelated with slower particles. This is a 
weaker condition than the "stosszahlansatz" that the two 
particle density be equal to a product of one-particle den- 
sities. 

Eq. reflects that large velocities undergo the fol- 
lowing cascade process 



V -> (pv, qv), 



(12) 



with the rate \v\ x . These cascade dynamics follow di- 
rectly from the collision rule Q by setting one of the 
incoming velocities to zero. Even though the number 
of particles is conserved, the number of energetic parti- 
cles doubles in each cascade event (Fig. Moreover, 
momentum is conserved but energy is dissipated in each 
cascade event: it is transferred from large velocities to 
smaller velocities. 



In general dimensions, the collision rule is 
vi = ui - (1 -p)(ui - u 2 ) • nn. 



(13) 



Here n = n/n with n = |n| is a unit vector parallel to 
the impact direction n (connecting the particle centers), 
Vi j2 are the post-collision velocities, and Ui j2 are the pre- 
collision velocities. The normal (to n) component of the 
relative velocity is reduced by the restitution coefficient 
r = 1 — 2p as follows, (vi — V2) • n = — r (ui — u 2 ) • n 
and the energy dissipated equals p(l — p)\(ui — u 2 ) • n| 2 . 
Similarly, the general collision rate J2J becomes 



if(ui,u 2 ) = |(ui 

/d(v) satisfies 



U2)-n| A . The velocity distribution 



~fc(v) 



dndui du 2 |(ui - u 2 ) • n| A / d (ui)/ d (u 2 ) 
x [<5(v-vi)-*(v-ui)]. (14) 



In addition to integration over the incoming velocities, 
an additional integration over the impact direction is re- 
quired, and this integration is normalized, J dh = 1. The 
impact angle is assumed to be uniformly distributed. 

The dynamics of large velocities v — -> 00 are simplified 
as in the one-dimensional case. The integration over the 
incoming velocities is separated into an integral over a 
small velocity and an integral over a large velocity u. 
The former integration is immediate, 

^/d(v) = J J dhdu\u ■ n| A / d (u) x (15) 

[<J(v— (1— p)u ■ nn) + i(v-u+(l-p)u ■ nn) + 5(v — u)] . 

Let /1 = (u • n) 2 ; in other words, if 9 is the an- 
gle between the dominant velocity and the impact 
angle, then /i = cos 2 8. There are two gain terms 
corresponding to the two cases v=(l-j))u-nn and 
v = u — (1 — p)u ■ nn. These collision rules, together 
with the impact angle, dictate the magnitudes of the 
post-collision velocity in terms of the pre-collision veloc- 
ity, v — au and v — f3u with the following stretching 
parameters 



a={l-p)^\ 
(3=[l-(l-p^] 1/2 . 



(16a) 
(16b) 



The parameter a follows from u = n and 
the parameter (3 is obtained by introducing 
w = v — u and then employing the collision rule 
w = — w • n = (1 — p)u • n = (1 — p)u/i 1 / 2 and the iden- 
tity v 2 = u 2 + w 2 — 2ww/i 1 ^ 2 . The integration over the 
large velocity u includes separate integrations over the 
velocity magnitude u and over the velocity direction 
u but since this angle is a unique function of the 
impact angle, the latter integration is immediate. Using 
the isotropic velocity distribution /<j(v) = SdV^ 1 f(v) 
with Sd J = 1 and Sd the area of the 

c?-dimcnsional unit hypersphere, Eq. (|15|l simplifies to 



-x dfjv) 
dt 



diidu\uiJ 1 ' 2 \ x u d - 1 f(u) 



,.V2|A, ( d-lfC,^ ( 17 ) 

x [S(v — au) + S(v — (3u) — S(v — u)] . 



Finally, we simplify the angular integration, 
dn oc sin d_2 6 d9. Denoting the angular integration 
with angular brackets (g) = Jrfn g(n), we have 




FIG. 2: The exponent a versus the restitution coefficient r for 
hard-spheres (A = 1) and Maxwell molecules (A = 0). The 
top two curves are for d — 3 and the bottom two curves are 
for d = 2. 



(g) = C / dMsC/z)/*- 1 /^!-^ 



(18) 



The constant C = l/B(\,^-), with B(a, b) the beta 
function, is set by normalization. The linear equation 
governing the tail of the distribution is therefore 



df(v) 
dt 

(19) 

As in one-dimension, large velocities undergo the cas- 
cade process 



v — > ( av , 0v ), 



(20) 



but in general dimension, the stretching parameters ac- 
quire a dependence on the impact angle. In each colli- 
sion, the total velocity magnitude increases, despite the 
fact that the total energy decreases, as reflected by the 
following two inequalities 



a + (3 > 1, 
a 2 +(3 2 < 1. 



(21a) 
(21b) 



Equalities occur in the limiting cases: the total velocity 
magnitude is conserved in one dimension where collisions 
are always head-on (jj, = 1) and, of course, the total en- 
ergy is conserved for elastic collisions (p = 0). Actu- 
ally, a stronger statement than (|21|) holds: the quantity 
M s (a, (3) = a s + (3 s — 1 is positive for s < 1, negative for 
s > 2, and it may be either positive or negative in the 
range 1 < s < 2, depending on the impact angle. 



V. STATIONARY STATES 

We have seen that the velocity distribution decays al- 
gebraically for one-dimensional Maxwell molecules. The 



linear equation for the tail of the distribution shows that 
this behavior extends to all A and all p in one-dimension. 
The power-law velocity distribution 



f(v) 



(22) 



satisfies the linear Boltzmann equation i|ll|) with the time 
derivative set to zero when p CT_1_A + g cr ~ A ~ 1 = 1. Since 
the collision parameters satisfy p + q = 1 , the character- 
istic exponent in one-dimension is simply 



a = 2 + X. 



(23) 



Of course, the power-law behavior applies only for the 
tail of the distribution. 

Algebraic behavior holds in arbitrary dimension. Sub- 
stituting Eq. 11' I'D into the general linear Boltzmann equa- 
tion i|19|) , the characteristic exponent is root of the equa- 
tion 



j-d-X 



(3° 



-d-X 



-l)/// 2 ) = 



(24) 



This transcendental equation can be re- written explicitly 
in terms of the gamma function and the hypergeometric 
function 441 



j+A- 



A+l 
2 ' 



d+X 



,i-p 2 ) _r(g=f±i)r(4±*) 



(i-pY 



-d-X 



r(§)r(Mi) 



(25) 



The exponent a = cr(d, A, r) varies continuously with the 
spatial dimension d, the homogeneity index A, and the 
restitution coefficient r (figure |2J). 

According to the bounds (12111 the left hand side of 
Eq. I|24|) is positive when a — d — A < 1 but negative 
when a — d — A > 2. Therefore, this quantity changes 
sign when 1 < u — d — A<2 leading to the relatively 
tight bounds 



l + A<er<d + 2 + A. 



(26) 



The lower bound l|23|) is realized in one-dimension where 
the collisions are always head-on, while the upper bound 
is approached, a — > d + 2 + A, in the quasi-elastic limit 
r — * 1. We note that the two limiting cases of one- 
dimension and elastic collisions do not commute. More- 
over, the zero dissipation limit is singular: Maxwellian 
distributions occur when the collisions are elastic [2jJ . 

Since the energy lost in each collision is proportional 
to (Ad) 2 and the collision rate is proportional to |Au| A , 
the energy dissipation rate is related to the following in- 
tegral, T ~ (v 2+x ) where (g(v)) = S d J dv v^ 1 f(v)g(v). 
Hence, the bound a < d + 2 + A implies that the total 
dissipation rate is divergent. This is a generic feature of 
the stationary solutions, and in fact it shows why Haft's 
cooling law dT/dt = — r, where T — (v 2 ) is the granular 
temperature, does not apply: this rate equation assumes 
finite dissipation rates. In contrast, the total energy may 
be either finite or infinite because both a > d + 2 and 
a < d+2 are possible. The stationary states studied here 
appear to be fundamentally different than the infinite en- 
ergy solutions of the elastic Boltzmann equation because 
they require dissipation and because they always involve 
infinite dissipation |45j. 

The characteristic exponent increases monotonically 
with the spatial dimension, the homogeneity index, and 
the restitution coefficient. Thus, fixing d and A, the com- 
pletely inelastic case (r — 0) provides a lower bound 
for a with respect to r (figure^. For hard-spheres the 
completely inelastic limit yields a = 4.1922 and a = 
5.23365 in two- and three-dimensions, while for Maxwell 
molecules the corresponding values are a — 3.19520 and 
a = 4.28807. 

The power-law behavior is in sharp contrast with the 
stretched exponential tails f(v) ~ exp(— \v\ 5 ) that typi- 
cally characterize granular gases. For freely cooling gases, 
S = A as in ®, and for thermally forced gases, S = l+A/2 
HH l48l |49| . Both behaviors immediately follow from 
the linear Boltzmann equation (|19fl : in the forced case, 
the time derivative in (|19fl is replaced by the diffusive 
forcing term d/dt — > Z?V 2 . Only in the limiting case of 
freely cooling Maxwell molecules do power-law velocity 
distributions arise, but the solutions are not stationary 
and the characteristic exponent differs from the station- 
ary solutions [5(J, |5ll |52, |53| • 



VI. DRIVEN STEADY-STATES 

In this section we describe driven, non-equilibrium 
steady-states that are identical, except for a high veloc- 
ity cut-off, to the stationary states described above. In 
these steady-states energy is injected at a large veloc- 
ity scale, cascades to small velocities, and is dissipated 
over a broad power-law range. The energy injection scale 
V must be well separated from the typical velocity scale 
i>o, but otherwise, the injection mechanism is not unique. 
We study several concrete cases where energy is injected 
with a small rate at a large velocity. 
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FIG. 3: The velocity distribution for one-dimensional 
Maxwell molecules. The solid line is a Lorentzian and the 
typical velocity is Vo = 0.055. 

As we have seen in the exactly soluble case of one- 
dimensional Maxwell molecules, there is a family of 
steady-state solutions characterized by the typical veloc- 
ity vq'- if f(v) is a steady-state solution, so is Vq f (v /vq) 
for arbitrary vq. Let the energy injection rate (per parti- 
cle) be 7 and let the injection velocity be V. This scale 
sets an upper cutoff on the velocity distribution, beyond 
which the distribution should rapidly vanish. Since the 
system is at a steady-state, the dissipation rate 

T ~ (v 2+x )~ J V dvv d+1+x v d f(v/v ) (27) 
~ V x+2 {V/v a ) d -°, 

must be balanced by the energy injection rate "fV 2 , lead- 
ing to a general relation between the injection rate, the 
injection velocity and the typical velocity, 

~f~V x (V/v ) d -°. (28) 

To verify the theoretical predictions, we performed 
Monte Carlo simulations. Collisions are simulated by 
selecting two particles at random with a probability pro- 
portional to the collision rate and then updating their 
velocities according to the collision rule <|13(l . Energy is 
injected with a small rate using the following "lottery" 
implementation. An energy loss counter keeps track of 
the cumulative energy loss. With a small rate, a ran- 
domly chosen particle is "awarded" an energy equal to 
the reading on the loss counter. Subsequently, the loss 
counter is reset to zero. With this protocol, the kinetic 
energy remains practically constant, and moreover, en- 
ergy injection occurs only at large velocity scales. For 
one-dimensional hard-spheres, we tested a different injec- 
tion mechanism. The injection energy was drawn from 
a Maxwell-Boltzmann distribution with a very large en- 
ergy. With a small rate, this energy was added to a 
randomly chosen particle. 

We simulated completely inelastic Maxwell molecules 
and hard spheres in one- and two-dimensions starting 




FIG. 4: The velocity distribution for Maxwell molecules. 
The top curves correspond to one-dimension and the bottom 
curves to two-dimensions. 

with a uniform velocity distribution with support in the 
range [—1 : 1]. After a short transient, the system reaches 
a steady-state. For the special case of one-dimensional 
Maxwell molecules, we verified that the velocity distri- 
bution is Lorentzian (figure 0J). In all cases, the tail of 
the velocity distribution decays as a power-law, and the 
exponent a is in excellent agreement with the theoretical 
prediction, Eq. I|25[) . Maxwell molecule simulation re- 
sults are displayed in figure 0] and hard sphere simulation 
results in figure 

The energy balance relation (|28|) . combined with the 
constant energy condition (v 2 ) ~ 1, imposed in our simu- 
lations, yields an estimate for the typical velocity. Differ- 
ent behaviors emerge for finite energy and infinite energy 
distributions. 

When a < d+2, the constant energy constraint implies 
yd+2-0 ^ v d-a ^ combined with energy balance 
reveals how the maximal velocity and the typical velocity 
scale with the injection rate 

y^ 7 -5 J ^ j (29a) 

v ~ 7 <--rfK2-A) _ (29b) 

Simulations with d = 1, A = 0, and 7 = 10~ 4 , are charac- 
terized by V « 10 2 and v « 10~ 2 , consistent with these 
scaling laws. 

In the complementary case, a > d + 2, the typical 
velocity vq ~ 1 is set by the initial conditions because 
{v 2 ) ~ Vq. Energy balance f2*%|l yields 

V^ 7 -?=h=x. (30) 

Simulations with d = 2, A = 1, and 7 = 10 -2 should be 
characterized by the injection scale V ~ 50, as in this 
case a = 4.15. The data is consistent with this estimate. 

As long as the system is sufficiently large, there is no 
dependence on the system size (the number of particles). 
The total energy and the total dissipation rate are pro- 
portional to the system size, and in general, all thermo- 
dynamic properties are extensive. 
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FIG. 5: The velocity distribution for hard spheres in 
one-dimension (top curves) and in two-dimensions (bottom 
curves) . 



Based on the theoretical and the simulation results, 
we conclude that there may be qualitative differences be- 
tween the finite energy and infinite energy cases, but that 
fundamentally, the steady-state solutions are of one na- 
ture. They represent a nonequilibrium stationary-state 
where energy is injected at velocity scale V and dissi- 
pated at velocity scale vq. These scales are set by the 
injection rate and the injection protocol. 



VII. TIME-DEPENDENT STATES 

What happens to these steady states when energy in- 
jection is turned-off? Steady-state solutions of the type 
(|22|l can be realized only up to some upper cutoff. Such 
truncated power-law distributions are still compact, and 
thus, in the absence of energy input, they should undergo 
free cooling with all energy eventually dissipated from the 
system. 

Therefore, we anticipate that there is a time-dependent 
velocity cut-off V(t). Below this scale the distribution is 
nearly the same as the stationary distribution but above 
this scale, the distribution has a sharp tail, analogous 
to the freely cooling state Q. Thus, the distribution is 
of the form f(v,t) = f(v,V) such that for v < V(t), 
f(v;V) « fs(v) while for v > V(t), the distribution de- 
cays faster than a power law. Here f s (v) is the stationary 
solution of the full Boltzmann equation. We assume that 
the cut-off scale is much larger than the typical veloc- 
ity, V > wo and, without loss of generality, set vo = 1. 
The assumption that the distribution is unmodified be- 
low the cut-off velocity is consistent with the character 
of the energy cascade. Furthermore, we expect that the 
functional form of the cut-off depends only on the scaled 
variable, v/V. 

First, consider the time dependence of the cut-off scale 
V(t). Given the assumption that cooling occurs only 
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through a decrease in the cut-off scale, the rate of change 
of the energy is 

dE d{v 2 ) d f VW ,_, d+1 



dt dt dt 

„ v d+1 -°— 

dt ' 



dvv a+L f s (v) (31) 



The decrease in energy equals the dissipation rate T ~ 
yd+2+x-a from Eq ip^ji^ showing that the cut-off veloc- 
ity obeys Haff 's cooling law, 



dV 
dt 



= -cV 



1+A 



(32) 



Therefore, the cut-off velocity decays with time as follows 



V(t) 



V x (0) 



l + cXV x (Q)t 



1/A 



(33) 



where V(0) is the initial value of the cut-off. 

Restricting our attention to one-dimension we seek 
similarity solutions of the type 



(34) 



Here, f s (v) is the stationary solution of Eq. J3J that de- 
cays as a power-law at large velocities f s (v) ~ Av~ 2 ~ x . 
The cut-off function approaches unity at small argu- 
ments, 4>{x) — > 1 as x — > so that the stationary solution 
is recovered for v -C V. 

Substituting the time dependent form (|34l) into the 
linear governing equation (|ll|l yields 



dV/dt 
' v 1+x 



4f(x) 



P0 



(35) 

Assuming the cut-off velocity satisfies Eq. I|32fl with con- 
stant of proportionality c, the scaling function satisfies 
the linear and non-local differential equation 



c<j>'(x) 



.A— 1 



«K*) 



(36) 



with the boundary conditions 4>(0) = 1 and 4>{x) — > as 
x — ► oo. Note that <fi(x) must be non-analytic at x — 
because all its derivatives vanish at x = since p + q = 1 
and 0(0) = 1. 

For large arguments, the last term on the right hand 
side dominates, and therefore, the tail of the distribution 
is a stretched exponential as in © 



exp {—Cx > 



(37) 



with C — (Ac) -1 for A > 0. In the limiting case A = 
all terms on the left-hand side are comparable and the 
tail is algebraic: 4>(x) ~ x~ a with cer = 1 — p a+1 — q a+1 . 
Thus, both the decay of the cut-off velocity and the 
tail behavior are as for ordinary freely cooling solutions, 



Eq. There is, however, a difference since the distri- 
butions considered here have two characteristic velocities, 
V and vo and it is only the upper cut-off, V that evolves 
in time. After V and vo become comparable, the behav- 
ior crosses over to the homogeneous cooling state 
with a single characteristic velocity, vq. 

We now focus on completely inelastic hard-spheres 
(A = 1 and p = q = 1/2) for which an exact solution 
is possible. Integrating Eq. (|3"6*jl and imposing (j>(0) = 1 
gives c = (1 — p 2 — q 2 ) / °° dx 4>(x), but since the cut-off 
scale V is defined up to a constant, we may set the in- 
tegral value, J °° dx <j>(x) = 1 leading to c = 1 — p 2 — q 2 . 
When p = q = 1/2 then c = 1/2 and Eq. (|36|l becomes 



(j)'(x) =2{4>{2x)-4>(x)]. 



(38) 



This equation can be solved using the Laplace transform 
h(s) — J dx e~ sx 4>(x), that satisfies the non-local equa- 
tion 



(2 + s)h(s) = 1 + h{s/2) 



(39) 



with the boundary condition h(0) = 1 set by the 
normalization. Since <j>{x) — > 1 as x — > then 
h(s) — > s -1 as s — > oo, so we make the transforma- 
tion h(s) = s- 1 [1 - g(s)} with g(0) = 1 and g'(0) = -1. 
The auxiliary function g(s) satisfies a recursion-like equa- 
tion g(s) = (l + s/2) 1 g(s/2). Solving iteratively and 
invoking <?(0) = 1, the solution is the infinite product 
g(s) = n^Li(l + s /2 n ) _1 , and the Laplace transform is 



h(s) = - s [l- || (40. 



Since the infinite product has a series of simple poles at 
s = —2~ n for every integer n > 1, the scaling function is 
a sum of exponentials 



4>(x) = a n exp (— 2 n x) 



(41a) 



(41b) 



fe=i 



with the coefficients obtained as the residues to the poles 
a n = lim s ^_2" [(s + 2")/i(s)]. In contrast with freely 
cooling states (@J, the scaling function <t>(x) can be ob- 
tained exactly. 

The Laplace transform conveniently yields the limiting 
behaviors of the scaling function. The simple pole closest 
to the origin reflects the tail behavior 



(x) ~ ai exp(— 2 x) 



(42) 



as x — > oo with a\ = 3.46275 obtained from Eq. I|41bj) . 
More interesting is the small x behavior, reflected by the 
large s behavior 



dx[l 



= s^gis) s^exp [-C(lns) 2 ] 




> 10 




FIG. 6: The scaling function <j>(x) versus x for A = 1/2 
(dashed line) and A = 1 (solid line). 



as s — > oo with C = (21n2) _1 . The function g(s) was 
estimated by replacing the infinite product with a finite 
product 



x i n * on 



(43) 



with n* = \xi2 s. Inverting the log- normal Laplace trans- 
form using the steepest descent method, the leading cor- 
rection to the scaling function is log-normal as well 



1 — 4>{x) ~ exp 



I ! lni 

x 



(44) 



as x -> with A = C/4 = (81n2) _1 . Thus, the scaling 
function is perfectly flat near the origin as all its deriva- 
tives vanish at x = (figure 0). Physically, the small x 
behavior shows that there is a sizable range of velocities 
for which the time-dependent velocity distribution 1341) 
coincides with the steady-state solution, f(v,t) ~ f 8 (v). 

The series solution (1411) can be straightforwardly gen- 
eralized to all A > 



${x) = J2 a n exp [-(2"x) A ], 
n=l 

oo _^ 

n y~~2m 



k=l 

k^n 



)X(n—k) ' 



(45a) 
(45b) 



Making the transformation y = x x and setting 
the proportionality constant c = A _1 2~ A such that 
c = (1 - 2- x ) J™dy(j)(y), Eq. jSHJ is generalized, 
<t>'{y) — 2 A [0(2 A y) — 4>(y)] . Consequently, the Laplace 
transform is obtained from Eq. 1)41) I) by replacing 2™ with 
2 A ", and repeating the steps leading to (|4*T|) gives (@SJ). 
Figure H3 shows the scaling function for A = 1/2. As A de- 
creases the cut-off becomes less sharp and the flat region 
near x = less broad. 



FIG. 7: The velocity distribution f(v, i) versus v. Shown is 
the steady-state distribution before the injection is turned- 
off (solid line) and at three consecutive and equally-spaced 
later times (circles, squares, diamonds) during free cooling. 
Also shown for reference is a dashed line with slope -3. The 
velocity is in arbitrary units. 



In summary, we find that there are time-dependent 
states associated with the stationary states. In these 
transient states, the velocity distribution is characterized 
by a cut-off velocity scale that decays with time accord- 
ing to Haff 's law. Below this velocity, the energy cascade 
is unaffected and the velocity distribution agrees with 
the stationary distribution but above this scale, the dis- 
tribution is exponentially suppressed. We relied only on 
the linear Boltzmann equation to derive a scaling form 
for the cut-off function. Of course, the full nonlinear 
equation J2J is still relevant as it governs the dynam- 
ics of small velocities via the stationary solution f s (v). 
This guarantees that the velocity distribution is properly 
normalized, and specifically, that the integral over small 
velocities remains finite. 

We numerically integrated the hard-sphere Boltzmann 
equation in one-dimension to verify these predictions. 
Velocity bins are kept, each with a double precision num- 
ber representing the number of particles within that ve- 
locity range. In the simulation, two velocity bins are 
chosen randomly with a rate proportional to the collision 
rate. When two bins "collide" , particles are transferred 
from each into target bins, determined by the collision 
rule (|TJ). 

We generated the stationary distribution by injecting 
energy at a fixed rate. This was done by uniformly re- 
moving particles from the distribution and re-injecting 
them according to a Gaussian distribution with a large 
characteristic velocity. Once the system reaches the sta- 
tionary state, we turn off the energy injection and observe 
the distribution f(v,t) as it cools. 

Figure0shows the driven steady-state distribution and 
the freely cooling distributions at three later times. The 
results verify that the steady-state has a power-law tail, 



/.(«) 



and that the freely cooling distributions 



are close to the steady-state distribution for sufficiently 




FIG. 8: The scaling function underlying the velocity distribu- 
tion. The velocity distributions in figure [7| were normalized 
by the stationary distribution as in Eq. ||3]J. The solid line 
is the theoretical scaling function l|410 . 

small velocities. Figure |S] shows the same three time- 
dependent distributions divided by the steady distribu- 
tion as in Eq. 1|34|) and rescaled by the cut-off velocity 
V(t) to collapse the data onto the theoretical prediction 

The time dependence of the cut-off velocity, given by 
Eq. (j22J, holds until V is order Vq = 1. Thus, the life- 
time of the collapsing power-law solution approaches a 
constant of order unity as V(0) becomes infinite. Dur- 
ing most of the time that the power-law is collapsing, 
V decays algebraically with time, V(t) ~ t~ x / x . Figure 
[5] shows the cut-off velocity versus time together with a 
fit to the form (|33|l with A — 1. We also checked that 
the tail of the cooling distribution is exponential. Wc 
conclude that for completely inelastic hard-spheres, the 
simulation results are in excellent agreement with the 
theoretical predictions. 



VIII. CONCLUSIONS 

In summary, we find a new class of steady-state and 
time-dependent states for inelastic gases. In the nonequi- 
librium steady-states, energy is injected at large veloci- 
ties, it cascades down to small velocities, and it is dissi- 
pated over a power-law range. Generically, the steady- 
state distributions have a power-law high-energy tail. 
The characteristic exponents were obtained analytically 
and they vary with the spatial dimension and the colli- 
sion rules. Formally, the stationary solutions are charac- 
terized by an infinite dissipation rate, while the energy 
density may be either finite or infinite. In an actual par- 
ticle system, these steady-states may be realized only up 
to the energy injection scale, so that all thermodynamic 
characteristics including the dissipation rate and the en- 
ergy density are finite. 

When injection is turned-off, the velocity distribution 
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FIG. 9: The cut-off velocity V(t) as a function of time t (cir- 
cles). The solid line is a fit to Eq. . Also shown is a 
broken line with slope —1. 



is stationary only in a shrinking range of velocities and 
it decays sharply as a stretched exponential at large ve- 
locities. For completely inelastic collisions, the scaling 
function underlying this behavior can be obtained exactly 
from the linearized Boltzmann equation and at small ve- 
locities, there is a subtle log-normal correction to the 
power-law behavior. Although we analyzed only the one- 
dimensional case, we expect the same behavior in higher- 
dimension. These time-dependent states can be loosely 
thought of as a hybrid between a steady-state solution 
and the well-known, freely-cooling solution. Both the 
time-dependence of the characteristic velocity and the 
decay at large velocities are similar, but not identical, to 
those occurring for freely cooling granular gases. After 
the cut-off velocity becomes comparable to the typical ve- 
locity, the velocity distribution presumably crosses over 
to the freely cooling solutions. 

Cascade processes occur in many other physical sys- 
tems. Mathematically, the inelastic cascade process for 
one-dimensional hard-spheres is identical to that found 
for the grinding process in Ref. |5J|, and in both prob- 
lems cr = 3. Indeed, the cascade process (I12fl is equiva- 
lent to a fragmentation process. In fluid turbulence, the 
fluid is forced at a large spatial scale, energy cascades 
from large scales to small scales, where it is dissipated 
due to viscosity [5f|- Actually, the situation found here 
for granular gases is analogous to wave turbulence, that 
is described by a kinetic theory for wave collisions |56| . 
One difference with the Kolmogorov spectra of fluid tur- 
bulence is that the characteristic exponents are irrational 
because they do not follow from dimensional analysis. 

Inelastic cascades are a direct consequence of the col- 
lision rule and they are described by a linear equation. 
This equation should be valid under very broad condi- 
tions and it can be generalized to nonuniform distribu- 
tions of impact angles and collision parameters as well as 
polydisperse granular media. Additionally, the cascade 
dynamics should extend to viscoelastic collision rules be- 
cause the restitution coefficient depends only weakly on 
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the relative velocity for energetic collisions [l|. 

The most significant condition for these steady-states 
concerns the driving mechanism: the injection rate must 
be small compared to the collision rate and the injec- 
tion energy large compared to the typical energy. We 
propose that stationary distributions may be achieved in 
driven granular gas experiments where energy is injected 
at very large velocity scales. Algebraic tails with expo- 
nents comparable with these reported here were observed 
recently in sheared granular layers, but these experiments 
involved frictional, rather than collisional, dynamics |57| . 

Inelastic cascades should also arise when an energetic 
particle hits a static medium of inelastic particles, or al- 
ternatively, a background of slowly moving particles. In- 
deed, the collision dynamics in this case reduce to the 
inelastic cascade discussed in this paper. This setup may 
be interesting to study theoretically and experimentally. 



In closing, the kinetic theory of inelastic collisions is 
remarkable as the nonlinear Boltzmann equation admits 
a number of distinct solutions including steady-states, 
transient states, and hybrid states that interpolate be- 
tween the two. Nonlincarity, nonlocality, and the lack of 
energy conservation are responsible for this remarkable 
complexity. We end with an open question: do other 
families of solutions exist? 
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